3d surface related multiple elimination for wide azimuth seismic data

ABSTRACT

According to a preferred aspect of the instant invention, there is provided herein a system and method for the imaging and monitoring of the subsurface hydrocarbon reservoirs and other subsurface features, wherein the collected seismic data contain multiples therein. In brief, the instant invention is a method of removing multiples using SRME, wherein wide azimuth seismic data are used to directly compute a multiple prediction via a mixed space-wavenumber-frequency domain implementation.

RELATED CASES

This application claims the benefit of U.S. Provisional PatentApplication Ser. No. 60/956,278 filed on Aug. 16, 2007, and incorporatessaid provisional application by reference into this disclosure as iffully set out at this point.

TECHNICAL FIELD

This invention relates to the general subject of geophysical explorationfor hydrocarbons and, in particular, to methods for collecting andanalyzing seismic data that evidence multiples therein.

BACKGROUND OF THE INVENTION

A seismic survey represents an attempt to image or map the subsurface ofthe earth by sending energy down into the ground and recording the“echoes” that return from the rock layers below. The source of thedown-going sound energy might come, for example, from explosions orseismic vibrators on land, or air guns in marine environments. During aseismic survey, the energy source is placed at various locations nearthe surface of the earth above a geologic structure of interest. Eachtime the source is activated, it generates a seismic signal that travelsdownward through the earth, is reflected, and, upon its return, isrecorded at a great many locations on the surface. Multiplesource/recording combinations are then combined to create a nearcontinuous profile of the subsurface that can extend for many miles. Ina two-dimensional (2D) seismic survey, the recording locations aregenerally laid out along a single line, whereas in a three dimensional(3D) survey the recording locations are distributed across the surfacein a grid pattern. In simplest terms, a 2D seismic line can be thoughtof as giving a cross sectional picture (vertical slice) of the earthlayers as they exist directly beneath the recording locations. A 3Dsurvey produces a data “cube” or volume that is, at least conceptually,a 3D picture of the subsurface that lies beneath the survey area. Inreality, though, both 2D and 3D surveys interrogate some volume of earthlying beneath the area covered by the survey.

A conventional seismic survey is composed of a very large number ofindividual seismic recordings or traces. In a typical 2D survey, therewill usually be several tens of thousands of traces, whereas in a 3Dsurvey the number of individual traces may run into the multiplemillions of traces. Chapter 1, pages 9-89, of Seismic Data Processing byOzdogan Yilmaz, Society of Exploration Geophysicists, 1987, containsgeneral information relating to conventional 2D processing and thatdisclosure is incorporated herein by reference. General backgroundinformation pertaining to 3D data acquisition and processing may befound in Chapter 6, pages 384-427, of Yilmaz, the disclosure of which isalso incorporated herein by reference.

A seismic trace is a digital recording of the acoustic energy reflectingfrom inhomogeneities or discontinuities in the subsurface, a partialreflection occurring each time there is a change in the elasticproperties of the subsurface materials. The digital samples are usuallyacquired at 0.002 second (2 millisecond or “ms”) intervals, although 4millisecond and 1 millisecond sampling intervals are also common. Eachdiscrete sample in a conventional digital seismic trace is associatedwith a travel time, and in the case of reflected energy, a two-waytravel time from the source to the reflector and back to the surfaceagain, assuming, of course, that the source and receiver are bothlocated on the surface. Further, the surface location of every trace ina seismic survey is carefully tracked. This allows the seismicinformation contained within the traces to be later correlated withspecific surface and subsurface locations, thereby providing a means forposting and contouring seismic data—and attributes extractedtherefrom—on a map (i.e., “mapping”).

Many variations of the conventional source-receiver arrangement are usedin practice, e.g. VSP (vertical seismic profile) surveys, ocean bottomsurveys, etc. Further, the surface location of every trace in a seismicsurvey is carefully tracked and is generally made a part of the traceitself (e.g., as part of the trace header information). This allows theseismic information contained within the traces to be later correlatedwith specific surface and subsurface locations, thereby providing ameans for posting and contouring seismic data—and attributes extractedtherefrom—on a map (i.e., “mapping”).

For all of the subsurface information that might be acquired via seismicdata, this method is not without its problems. For example, oneparticularly troublesome problem in seismic data collection and analysisis the identification and removal of multiples. Those of ordinary skillin the art will understand that multiples in seismic data occur when therecorded seismic data contain energy that has been reflected more thanonce in the subsurface. Multiples often appear to all intents andpurposes to be valid seismic reflectors and, to the extent that they areinterpreted as such, can give rise to incorrect interpretations of thesubsurface layer configuration, thereby potentially resulting in thedrilling of dry holes.

One method of attenuating multiples that has had some success is knownas “surface-related multiple elimination” or “SRME”, hereinafter. Inbrief, this method operates by creating a dataset that contains onlypredictions of the multiples that are present in the data. Specifically,the method seeks to predict the seismic expression of multiples thathave experienced one or more reflections between, for example, theair-water interface and the subsurface. Then, the predicted multiplesare subtracted from the original data leaving behind (at leasttheoretically) only the primary energy.

The application of SRME has provided significant added value topractical processing projects over the last decade, yet the promise of atrue 3D solution has been elusive. It has long been known that full 3DSRME multiple prediction is best applied to data coverage with unaliasedsource and receiver sampling over a full range of azimuths and offsetswhich is only roughly approximated by the conventional narrow azimuthtowed streamer (NATS) 3D seismic data typically available today. In moreparticular, and has been observed in a variety of different contexts, inorder to accurately image complex subsurface structures the structureshould be illuminated from a variety of different offsets and azimuths.Wide azimuth surveys have been done for many years and have often provento yield superior data that can be subsequently migrated or otherwiseimaged to produce an improved picture of the subsurface as compared withtraditional/narrow azimuth surveys. Where the survey is not designed toacquire wide azimuth information, methodologies designed to estimate therequired wide azimuth seismic data from narrow azimuth data haveimproved the ability to predict complicated out-of-plane multiples(e.g., diffractions) to a limited degree.

Conventional methods for predicting and attenuating surface-relatedmultiples either assume that the subsurface is 2D and thus areapplicable only to 2D data shot in the dip direction or they assume thatwide azimuth seismic data can be adequately extrapolated from narrowazimuth data to enable a 3D multiple prediction.

Heretofore, as is well known in the geophysical prospecting andinterpretation arts, there has been a need for a method of using seismicdata to obtain image of the subsurface that does not suffer from thelimitations of the prior art. Accordingly, it should now be recognized,as was recognized by the present inventors, that there exists, and hasexisted for some time, a very real need for a method of geophysicalprospecting that would address and solve the above-described problems.

Before proceeding to a description of the present invention, however, itshould be noted and remembered that the description of the inventionwhich follows, together with the accompanying drawings, should not beconstrued as limiting the invention to the examples (or preferredembodiments) shown and described. This is so because those skilled inthe art to which the invention pertains will be able to devise otherforms of this invention within the ambit of the appended claims.

SUMMARY OF THE INVENTION

According to a preferred aspect of the instant invention, there isprovided herein a system and method for the imaging and monitoring ofthe subsurface hydrocarbon reservoirs and other subsurface features,wherein the collected seismic data contain multiples therein. In brief,the instant invention is a method of removing 3D multiples using SRME,wherein wide azimuth data seismic are used to directly compute the 3Dmultiple prediction via a mixed space-wavenumber-frequency domainimplementation; no direct data extrapolation is used as part of theprocess.

The instant invention teaches a method of predicting free-surfacemultiples for wide or multi-azimuth seismic data. The invention ispreferably implemented in the space-wavenumber-frequency domain. Thisprovides for the efficient computation of predictions for large datavolumes, output of the predicted multiples directly onto the acquisitionlocations, and sources and receivers that need not be explicitlyco-located at or near the surface.

Wide-azimuth towed streamer data presents an opportunity to improve theaccuracy of 3D SRME multiple prediction compared to methods that usenarrow-azimuth towed streamer data alone. A mixed domain formulationsuch as that taught herein is well suited for computing 3D multiplepredictions for wide azimuth data. Further, no one has heretoforeutilized a mixed space-wavenumber-frequency domain approach to 3D SRME.

The foregoing has outlined in broad terms the more important features ofthe invention disclosed herein so that the detailed description thatfollows may be more clearly understood, and so that the contribution ofthe instant inventor to the art may be better appreciated. The instantinvention is not to be limited in its application to the details of theconstruction and to the arrangements of the components set forth in thefollowing description or illustrated in the drawings. Rather, theinvention is capable of other embodiments and of being practiced andcarried out in various other ways not specifically enumerated herein.Finally, it should be understood that the phraseology and terminologyemployed herein are for the purpose of description and should not beregarded as limiting, unless the specification specifically so limitsthe invention.

BRIEF DESCRIPTION OF THE DRAWINGS

Other objects and advantages of the invention will become apparent uponreading the following detailed description and upon reference to thedrawings in which:

FIG. 1 illustrates the general environment of the instant invention

FIG. 2 contains a preferred logic that describes one aspect of theinstant invention.

FIGS. 3A and 3B illustrate schematically how Fourier coefficients from ahorizontally transformed shot and receiver record are located.

FIG. 4 contains a schematic illustration of a typical seismic processingprocess.

DETAILED DESCRIPTION

While this invention is susceptible of being embodied in many differentforms, there is shown in the drawings, and will herein be described,some specific embodiments of the instant invention. It should beunderstood, however, that the present disclosure is to be considered anexemplification of the principles of the invention and is not intendedto limit the invention to the specific embodiments or algorithms sodescribed.

GENERAL ENVIRONMENT OF THE INVENTION

FIG. 1 illustrates the general environment in which the instantinvention would typically be used. As a first step, a seismic surveywill be designed (step 110), wherein the survey geometry, sample rate,number of sources/receivers, etc. would typically be selected in orderto image a preferred subsurface target. Among the many parameters thatmight be considered in formulating the survey design are:

-   The surface area to be covered by the survey;-   Whether the survey will be conducted on land, offshore, or some    combination of the two environments;-   The depth of the target;-   The 3-D structure of the target (including its 2-D or 3-D dip, if    any);-   Whether the design will utilize an “end on” configuration (wherein    all of the active receivers are on the same side of the source) or a    “split spread” configuration (i.e., wherein active receivers are    placed both ahead of and behind of the source);-   The maximum offset (i.e., in the case where an active source is used    the distance from the source to the most distant active receiver)    and minimum offset (i.e., the distance from the source to the    closest active receiver);-   The receiver-to-receiver spacing;-   The source-point spacing if a controlled source is used (i.e., the    shot-to-shot spacing, where “shot” is used in the sense of “source    activation point”);-   The relation between source-points and receiver-points (e.g.,    sources near to receivers, sources midway between receivers, etc.);-   The frequencies expected in the received data;-   The strength of the sources, and the sensitivity of the receivers,    etc.    Of course, choosing survey parameters such as these for a seismic    survey is routinely done with an eye to obtaining the best possible    data and such design choices that are well within capability of one    of ordinary skill in the art. Further, those of ordinary skill in    the art will recognize that many of the previous parameters are    interrelated (e.g., specification of the target depth determines in    a general way a preferred maximum offset, etc.).

Next, equipment (including geophones and/or hydrophones or other seismicreceivers, as well as recording instruments, etc.) will be typicallymoved to and set up in the field at least approximately according to theplanned survey design 110. Next, and as is described more fully below, asurvey will be conducted 120 that is preferably at least approximatelyin accordance with the original design. Of course, it is certainlypossible that on-site changes will need to be made to the surveyas-designed. However, generally speaking the goal of the field crew isto replicate the parameter choices of the survey designer as closely asis possible.

After positioning the source and receivers, the data will preferably becollected conventionally depending on the sort of survey that is beingtaken. For example, if a controlled source survey is conducted eachsource activation might be accompanied by 8 seconds or so of recordingat a 2 ms sample interval, with the exact length of each recording andsample rate being depending on a number of factors well known to thoseof ordinary skill in the art.

As is typical in seismic surveys, the source will be activated and theresulting seismic signals sensed by the receivers and converted toelectrical energy which is subsequently digitized and recorded. Theresponse of each receiver or receiver group to the source willpreferably be captured digitally as a function of time and stored onoptical, magnetic, or other media in computer readable form fortransportation to a centralized computing facility where the data willbe processed, interpreted, and integrated into other data taken over thesame prospect. That being said, in some instances some amount of initialprocessing 130 will be applied to the data while it is in the field. Forexample, such in-field processing might be done in order to verify thequality of the data that are being collected. In other instances, thedata might be processed to see whether or not the target subsurface rockunits are being imaged adequately. In any case, after field processingthe data will usually at least be relatable to specific locations on thesurface of the earth.

Although the data that are collected according to the instant inventionmay be processed to some extent in the field (step 130), eventually itwill typically be transferred to a processing center where morecomputing resources 150 and algorithms 140 are available. In theprocessing center a variety of processes 140 might be applied to thedata to make them ready for use by the explorationist. At some point theprocessed data traces will likely be stored, by way of example only, onhard disk, magnetic tape, magneto-optical disk, DVD disk, or other massstorage means.

The methods disclosed herein would best be implemented in the form of acomputer program 140 that has been loaded onto a programmable computer150 where it is accessible by a seismic interpreter or processor. Notethat a programmable computer 150 would typically include, in addition toworkstations, mainframes, and super computers, networks of individualcomputers that provide for parallel and massively parallel computations,wherein the computational load is distributed between two or moreprocessors. As is also indicated in FIG. 1, in some preferredembodiments a digitized zone of interest model 160 would be specified bythe user and provided as input to the processing computer program. Inthe case of a 3-D geophysical data, the zone of interest 160 wouldtypically include specifics as to the lateral extent and thickness(which might be variable and could be measured in time, depth,frequency, etc.) of a subsurface interval. The exact means by which suchzones are created, picked, digitized, stored, and later read duringprogram execution is unimportant to the instant invention and thoseskilled in the art will recognize that this might be done any number ofways.

The algorithms that are used to process the seismic data might beconveyed into the computer 150 by means of, for example, a floppy disk,a magnetic disk, a magnetic tape, a magneto-optical disk, an opticaldisk, a CD-ROM, a DVD disk, a RAM card, flash RAM, a RAM card, a PROMchip, or loaded over a network.

Similarly, algorithms to implement a preferred embodiment of the instantinvention 145 will preferably be loaded into the computer 150 by anyconventional or unconventional means of conveying computer source orobject code.

After the seismic data has been subjected to the processes discussedherein, the resulting information will likely be displayed either on ahigh-resolution color computer monitor 170 or in hard-copy form as aprinted section or a map 180. The geophysical interpreter would then usethe displayed images to assist him or her in identifying subsurfacefeatures conducive to the generation, migration, or accumulation ofhydrocarbons.

Turning next to FIG. 4, as was indicated previously the instantinvention will preferably be made a part of and incorporated into aseismic processing sequence. Those of ordinary skill in the art willrecognize that the processing steps illustrated in FIG. 4 are onlybroadly representative of the sorts of processes that might be appliedto seismic data and the choice and order of the processing steps, andthe particular algorithms that have been selected, may vary markedlydepending on the individual seismic processor, the signal source(dynamite, vibrator, etc.), the survey location of the data (land, sea,etc.), the company that processes the data, etc.

As a first step, and as is generally illustrated in FIG. 4, a 2D or 3Dseismic survey is conducted over a particular volume of the earth'ssubsurface (step 410). The data collected in the field consist ofindividual unstacked (i.e., unsummed or prestack) seismic traces whichcontain digital signals representative of the volume of the earth thathas been surveyed. Methods by which such data are obtained and processedinto a form suitable for use by seismic processors and interpreters arewell known to those of ordinary skill in the art.

The purpose of a seismic survey is to acquire a collection of spatiallyrelated seismic traces over or proximate to a subsurface target of somepotential economic importance. Data that are suitable for analysis bythe methods disclosed herein might consist of, for purposes ofillustration only, an unstacked 2-D seismic line, an unstacked 2-Dseismic line extracted from a 3D seismic survey or, preferably, anunstacked 3D portion of a 3D seismic survey. The invention disclosedherein is most effective when applied to a group of unstacked seismictraces that have an underlying spatial relationship with respect to somesubsurface geological feature. Again, for purposes of illustration only,the discussion that follows will be couched in terms of traces containedwithin a 3-D survey, although any assembled group of spatially relatedseismic traces could conceivably be used.

After the seismic data are acquired, they are typically taken to aprocessing center where some initial or preparatory processing steps areapplied to them. As is illustrated in FIG. 4, a common early step 415 isdesigned to edit the input seismic data in preparation for subsequentprocessing (e.g., digitization, demux, wavelet shaping, bad traceremoval, etc.). This might be followed by specification of the geometryof the survey (step 420) and storing of a shot/receiver number and asurface location as part of each seismic trace header. Once the geometryhas been specified, it is customary to perform a velocity analysis andapply an NMO (normal move out) or similar correction to adjust eachtrace in time in order to account for signal arrival time delays causedby offset. Note, of course, that some or all of these steps (and thosethat follow) might be performed in the field before the data are broughtback to the processing center where greater computer resources aretypically available.

After the initial pre-stack processing has been completed, it iscustomary to condition the seismic signal that is located within theunstacked seismic traces before creating stacked (or summed) datavolumes (step 430). In FIG. 4 step 430 contains a typical “SignalProcessing/Conditioning/Imaging” processing sequence, but those skilledin the art will recognize that many alternative processes could be usedin place of the ones listed in the figure. In any case, the ultimategoal from the standpoint of the explorationist is generally theproduction of a stacked seismic volume or, in the case of 2D data, astacked seismic line for use in the exploration for hydrocarbons withinthe subsurface of the earth. As will become apparent hereinafter, themethods of the instant invention would most likely be applied as inconnection with step 430, i.e., multiple attenuation or removal.

As is suggested in FIG. 4, any digital sample within a stacked seismicvolume is uniquely identified by an (X, Y, TIME) triplet, with the X andY coordinates representing a position on the surface of the earth, andthe time coordinate measuring a recorded arrival time within the seismictrace (step 440). For purposes of specificity, it will be assumed thatthe X direction corresponds to the “in-line” direction, and the Ymeasurement corresponds to the “cross-line” direction, as the terms“in-line” and “cross-line” are generally understood in the art. Althoughtime is a preferred and most common vertical axis unit, those skilled inthe art understand that other units are certainly possible mightinclude, for example, depth or frequency. Additionally, it is well knownto those skilled in the art that it is possible to convert seismictraces from one axis unit (e.g., time) to another (e.g., depth) usingstandard mathematical conversion techniques.

The explorationist may do an initial interpretation 450 of the resultingstacked volume, wherein he or she locates and identifies the principalreflectors and faults wherever they occur in the data set. This might befollowed by additional data enhancement 460 and/or attribute generation(step 470) of the stacked or unstacked seismic data. In many cases theexplorationist will revisit his or her original interpretation in lightof the additional infonnation obtained from the data enhancement andattribute generation steps (step 480). As a final step, theexplorationist will typically use information gleaned from the seismicdata together with other sorts of data (magnetic surveys, gravitysurveys, satellite data, regional geological studies, well logs, wellcores, etc.) to locate subsurface structural or stratigraphic featuresconducive to the generation, accumulation, or migration of hydrocarbons(i.e., prospect generation 490). Of course, the ultimate goal is theproduction of hydrocarbons from the subsurface via a well or othermeans.

Preferred Embodiments

According to a first preferred embodiment, there is provided a systemand method for removing multiples in seismic data using SRME, whereinwide azimuth seismic data are used to directly compute the multiplepredictions via a mixed space-wavenumber-frequency domainimplementation. In the preferred embodiment, no direct dataextrapolation is used as part of the process.

The instant invention relates to the general subject of seismicexploration and, in particular, to methods for removing 3D multiplesassociated with the surface of the earth. The method is applicable to,but not limited to, wide-azimuth towed streamer seismic data.

Broadly speaking, according to the instant method, 3D multiples arepredicted and then subtracted from the measured data which contains theprimaries and multiples reflections. In the prior art multipleprediction is implemented either entirely in the space domain or thewavenumber domain. One innovative aspect of the instant invention is theimplementation of multiple prediction in a mixed space and wavenumberdomain. This has several advantages over typical implementations and/or3D multiple prediction methods that use a small subset of a wide azimuthdataset. Among these are:

-   1) The multiples are preferably predicted directly onto the acquired    data locations rather than onto an idealized grid;-   2) The wavenumbers at the surface can be limited to mitigate    aliasing effects;-   3) The 3D obliquity factor can be easily incorporated compared to    the space-only domain implementation;-   4) The sources and receivers do not need to be explicitly co-located    at or near the surface; and,-   5) The data need not be sorted to compute the multiple prediction.

In a preferred embodiment, the instant algorithm proceeds as issummarized below:

-   1) Transform the seismic data to the temporal frequency domain using    a time-direction Fourier transform;-   2) Read the entirety of the acquired seismic data for a single    temporal frequency;-   3) Select an inline and a crossline wavenumber;-   4) Separately transform the data over source and receiver domains    for the particular inline and crossline wavenumber;-   5) Select a particular output trace location (source and receiver    location);-   6) Compute the multiple prediction for that output location, the    inline and crossline wavenumber and the temporal frequency;-   7) Repeat steps 5 and 6 until there are no further output locations;-   8) Repeat steps 3-7 for a range of wavenumbers, summing the results;-   9) Output all of the output traces for a single frequency;-   10) Repeat steps 2-9; and,-   11) Inverse transform the output data back to the time domain.    Additional details related to the foregoing will be discussed below    in connection with FIG. 2.

In terms of equations, a preferred approach to the WATS 3D SRME problemsolved herein is to use a mixed wavenumber formulation to compute themultiple prediction. The conventional (prior art) equation forpredicting 3D multiples in the space-frequency domain can be written as:

${{M\left( {x_{g},\left. y_{g} \middle| x_{s} \right.,{y_{s};\omega}} \right)} = {\sum\limits_{y}\; {\sum\limits_{x}\; {{D\left( {x_{g},\left. y_{g} \middle| x \right.,{y;\omega}} \right)}{D\left( {x,\left. y \middle| x_{g} \right.,{y_{g};\omega}} \right)}}}}},$

where M is the multiple prediction data, D represents the input data, xand y refer to the trace coordinates, the subscripts s and g refer tosource and receiver respectively, and ω is the angular temporalfrequency. This equation clearly shows the need for source and receiverseverywhere on the acquisition surface. As was stated above, the currentWATS data does not lend itself to the straightforward application of theprevious equation.

The instant inventors have determined that the previous equation can bemodified to yield a mixed domain formulation, which reformulation is akey representation of the methods taught herein:

$\begin{matrix}{{{M\left( {x_{g},\left. y_{g} \middle| x_{s} \right.,{y_{s};w}} \right)} = {\sum\limits_{k_{x}}\; {\sum\limits_{k_{y}}\; {{D\left( {x_{g},\left. y_{g} \middle| k_{x} \right.,{k_{y};w}} \right)}2\; k_{z}{D\left( {k_{x},\left. k_{y} \middle| x_{s} \right.,{y_{s};w}} \right)}}}}},} & (1)\end{matrix}$

where k_(x) and k_(y) refer to inline and crossline wavenumbers andk_(z) represents the vertical wavenumber (w²/c²−k_(x) ²−k_(y) ²), wherec is the wave velocity at the near surface. The obliquity factor 2ik_(z)is often omitted in 2D and quasi-3D SRME formulations since theimprovement in accuracy obtained by including it can be smaller than thekinematic errors produced by inadequate data coverage. In addition, theobliquity factor is often more expensive to incorporate into thespace-frequency formulation. However, in the preferred embodiment, whenthe instant invention is utilized with WATS data, this factor willpreferably be included in the multiple estimation method taught herein,although it is not strictly necessary to do so.

Turning next to a detailed discussion of some preferred steps of theinstant invention, attention is drawn to FIG. 2. As an initial matter,note that the embodiment of FIG. 2 has been organized for purposes ofclarity in explaining the instant invention and not for purposes ofillustrating the most computationally efficient means of calculation.Those of ordinary skill in the art will readily understand how thecomputational scheme disclosed in this figure could be readily modifiedto improve its overall efficiency.

Turning to FIG. 2 in detail, as a first preferred step, a computerprogram will be initialized 200 according to methods well known to thoseof ordinary skill in the art (e.g., a program will be loaded intomemory, variables will be zeroed, files opened, etc.).

As a next preferred step 202, an unstacked seismic survey will beaccessed. In some preferred embodiments, the data might be storedlocally in memory (e.g., RAM) or on magnetic disk (or optical disk,tape, etc.). In other embodiments, the survey could be stored remotelyand accessed via a network. In any case, in the preferred arrangementseach seismic trace in the survey that is to be processed via the instantmethod will be stored as digital information in computer readable formaccording to a known format (e.g., the data will preferably be stored asindividual seismic traces comprising a trace header together withdigital signal values). Additionally, it is preferable that the surveybe a land or marine 3D survey.

According to a next preferred step 204, the seismic data will betransformed via Fourier transform in the time direction according tomethods well known to those of ordinary skill in the art. Thecalculation might be via a fast Fourier transform, a discrete Fouriertransform, etc.

Additionally, it should be noted and remembered that the Fouriertransform is a specific example of a class of data transformationsknown, for purposes of the instant disclosure, as projection transforms,i.e., transforms that calculate projections of data onto a (usually)infinite set of basis functions. More specifically, the Fouriertransform is an orthonormal projection transform because its basisfunctions are orthogonal and normalized to be of unit length/size. Otherexamples of projection transforms include the cosine transform, sinetransform, Walsh transform, wavelet transforms, the Radon transform,etc. Some of these (e.g., the Fourier, cosine, sine, and Walsh)transforms are usually chosen to be orthonormal, whereas others (e.g.,the wavelet and Radon transforms) may include basis functions that arenot normalized or even be mutually orthogonal. Thus, when the termprojection transform is used herein, that term should be interpreted tomean the Fourier transform (in the preferred embodiment) or, moregenerally, any projection transform of the sort discussed previously.Similarly, when the term projection transform coefficient is usedherein, that term should be interpreted to mean a Fourier transformcoefficient (in the preferred embodiment) or, more generally, anycoefficient obtained by projection (e.g., via least squares orotherwise) of the seismic data onto a particular basis function.

As a next preferred step, a shot record (or equivalently a shot gather)that is comprised of the Fourier transformed seismic traces will beaccessed and processed by calculating inline and crossline (2D) Fouriertransforms (step 206) separately for each frequency (horizontal) slice,a shot record typically being understood to be a collection of seismictraces that were recorded from the same shot or source excitation. SeeFIG. 3 and its discussion below for additional details of thiscomputation.

Note that the tranform directions “inline” and “crossline” are given asexamples only. Those of ordinary skill in the art will recognize thatany two (preferably) orthogonal directions could be used includingcalculating the 2D transform according to an arbitrary “x” and “y”coordinate system. Additionally, it should be noted that in order toperform this calculation is it not necessary for an entire shot recordto actually be read into memory at one time (i.e., an actual shot gatherneed not be formed by sorting the input data as is conventionally done)and the same statement applies to the receiver gathers discussed below.Methods of incrementally reading and processing such datasets so as tooperate in the shot, receiver, or other domain without sorting are wellknown to those of ordinary skill in the art. Thus, whenever it is saidherein that a shot, receiver, CMP, etc., record is read or accessed,that language should be broadly interpreted to include instances wherethe subject data traces are sorted and then read and held entirely inmemory as well as those instances where the data traces are arrangedaccording to some other ordering and are accessed sequentially or inparallel. Additionally, it should be noted that it may be necessary inthis or the next to create interpolated infill traces to reduce aliasingas is conventionally done.

The Fourier transform that is performed in connection with step 206 ispreferably applied to a horizontal slice through the unstacked shotvolume at each Fourier frequency, thereby transforming each slice intothe wavenumber (k_(x) and k_(y)) domain, i.e., a 2D transform of theconstant frequency slice will preferably be calculated. This processwill preferably be performed for every shot record in the survey. Noteonce again, as was explained previously, that in some preferredembodiments the individual shot records will not be assembled in memoryor on disk as there are numerically efficient means of calculating thequantities required herein without sorting the seismic dataset into shotorder. After preferably applying this computation to every shot recordin the survey, a transformed shot dataset (i.e., the “T_SHOT” dataset inFIG. 2) will be obtained.

As a next preferred step 208, Fourier transformed receiver records(i.e., all of the traces that were recorded by a single receiver) willbe accessed and a (horizontal or slice) Fourier transform calculated ateach frequency. As before, this process will be preferably repeated forevery receiver record in the survey and all frequencies, yielding theT_REC dataset of FIG. 2).

Preferably a loop will next be entered that selects output tracelocations (steps 210 through 232, inclusive). Note that the outputlocation 210 will preferably coincide with a seismic trace from thesurvey, but that is not strictly required.

Next, given the location of an output trace, that trace will beidentified within the T_SHOT and TREC datasets (step 212).

As a next preferred step, a loop over frequency will be entered (steps214 through 228, inclusive), wherein every Fourier frequency (or aselected range of frequencies) will be processed.

Within the frequency loop, the frequency slices from the T_SHOT andT_REC database will preferably be read (step 216) into memory, if theyhave are not already present there.

Next, preferably the instant invention will loop over inline andcrossline wave number (steps 218 through 224, inclusive). Note that thisloop is actually two dimensional in its preferred embodiment (i.e., adouble loop over k_(x) and k_(y) per equation (1) discussed previously)and those of ordinary skill in the art will readily be able to implementsuch. In some preferred embodiments, the range of each loop might gofrom negative to positive Nyquist according to a predeterminedwavenumber increment (which might potentially be different for theinline and crossline loops). As a specific example and for purposes ofillustration only, if the receiver cable were 2000 meters long withreceiver groups at intervals of 100 meters, the Nyquist wavenumber wouldbe 1/(2*100) or 1/200. In such a case, an increment approximately equalto 1/(length of the cable) or 1/2000 in this example would be apreferred initial estimate of this parameter value.

Note that there will likely be a tradeoff between the size of thewavenumber increment and the range of wavenumbers considered within thedouble loop and the quality of the end product. That is, smallincrements (and wavenumber ranges from 0 to Nyquist) will generallyyield at least somewhat better results as compared with largerincrements and a restricted range of wavenumbers. However, it is to beexpected that the computer run time will increase as the incrementdecreases. Similarly, if greater wavenumber ranges are used in thecalculation, the computer run time will also tend to increase. Those ofordinary skill in the art will be readily able to select theseparameters to obtain satisfactory results in a particular case and,indeed, in some instances these parameter choices will need to beselected on a trial-and-error basis.

Preferably, the instant invention will next identify a Fourier transformcoefficient at a given inline/crossline location in the T_SHOT dataset(step 218) and a corresponding Fourier coefficient at the same (orapproximately so) inline/crossline location within the T_REC dataset(step 220). Of course, if an exact match cannot be located because, forexample, of different shot and receiver spacing, such could readily beobtained by interpolation from nearby values.

FIGS. 3A and 3B contain a schematic representation of the previouscomputational steps. In the preferred embodiment, the selected outputtrace location 320 will be located within the transformed shot gather310 (T_SHOT) and receiver gather 315 (T_REC) calculated in steps 206 and208, respectively. At the current frequency w₀ in FIG. 3, 2D Fouriertransformed slices 340 and 345 will be identified and corresponding shot350 and receiver 355 transform coefficients (i.e., D(k_(x), k_(y)|x_(s),y_(s); w) and D(x_(g), y_(g)|k_(x), k_(y); w) in equation (1)) withineach of these slices will be identified and extracted.

Next, and preferably, the two Fourier transform coefficients 350 and 355that have been extracted from the shot 310 and receiver 315 gatherrespectively will be multiplied together, with the product preferablybeing multiplied by the obliquity factor 2ik_(z) according to equation(1). The product of these quantities will then preferably be added to anaccumulator of some sort (e.g., a memory location) per step 222 and anext inline/crossline wavenumber processed within the current frequencyslices (the “YES” branch of decision item 224).

After all of the wavenumbers have been processed and accumulated forthis frequency (the “NO” branch of decision item 228), preferably theaccumulated value for the current frequency will be stored for thecurrent trace (step 226) after which additional frequencies will besimilarly processed, with the accumulation trace, which contains one ormore frequencies therein, being written to storage (step 230) andidentified (e.g., by creating the appropriate trace header entry) withthe chosen output trace location per step 210. Of course, the outputtrace location could be in terms of inline/crossline trace numbers,latitude/longitude, or any arbitrary X/Y coordinate system and traceheaders are well suited to associate such infonnation with a seismictrace.

Obviously, the instant process could be repeated for any arbitrarynumber of output trace locations (e.g., the “YES” branch of decisionitem 232).

Finally, note that the output of the previous steps of the previousequation produces a frequency domain multiple prediction for eachselected location. Once one or more such multiple estimates have beencalculated, an inverse Fourier transform would typically be applied totransform the output traces from step 230 from the frequency back to thetime domain, at which point the time-domain multiple estimates might besubtracted from (or otherwise used to adjust) the original seismictraces, thereby producing seismic traces in which the multiples havebeen attenuated or eliminated. Additionally, the multiple estimateitself has independent values and might be studied in its own right todetermine the pattern of multiples in the data for purposes of trying todifferentiate between primaries and multiples in an unstacked or stackedseismic dataset.

Although the invention disclosed herein has been discussed in terms ofseismic traces organized (or organizable) into shot or receiver recordsor “gathers, that was done for purposes of specificity only and not outof any intent to limit the instant invention to operation on only thatsort of trace arrangement. So, within the text of this disclosure, theterms shot and receiver gather are used in the broadest possible senseof that term, and is meant to apply to seismic traces obtained fromconventional or unconventional 2-D and 3-D gathers, the most importantaspect of a “gather” being that it represents an organized collection ofunstacked seismic traces from either a 2-D or 3-D survey all of whichhave at least one shot location or receiver location at leastapproximately in common. This definition is specifically intended toencompass surveys where multiple shots are taken at the same location.

Additionally, in the previous discussion, the language has beenexpressed in terms of operations performed on conventional seismic data.But, it is understood by those skilled in the art that the inventionherein described could be applied advantageously in other subject matterareas, and used to locate other subsurface minerals besideshydrocarbons, e.g., coal. By way of additional examples, the sameapproach described herein could be used to process and/or analyzemulti-component seismic data, shear wave data, or model-based digitalsimulations of any of the foregoing.

It should also be noted that although the phrase “surface survey” hasbeen utilized herein to identify preferred embodiments, that term shouldbe broadly construed to include land surveys wherein the source andreceivers are located on the earth's surface, as well as marine surveyswhere the source/receivers are towed beneath the surface of the water,as well as marine surveys where the source is towed and the receiversare positioned on the bottom of the body of water (e.g., ocean bottomsurveys), as well as land surveys where the sources or receivers arepositioned at some predetermined depth, etc. Further, in some instances(e.g., VSP surveys) no surface survey will need to be conducted.

Further, although the preferred use of the instant invention is withactual seismic field data, it would also be suitable for work withmodel-generated data (e.g., synthetic data).

Additionally, it should be noted that while the instant inventionoperates best on seismic that was obtained from a wide azimuth survey,it would also be useful when used in conjunction with any sort ofseismic survey, whether or not collected with a view toward capturingwide azimuth information. Still, in the preferred embodiment a wideazimuth survey will be collected, as it will typically provide betterresults than one with a narrow range of azimuths. Methods of designingand acquiring wide azimuth surveys are well known to those of ordinaryskill in the art.

Still further, it should be noted and remembered that because of thesimple relationship between incoming/upcoming waves and downgoing/outgoing waves at the free surface, the instant invention could beimplemented, for example in the ω-k domain (i.e., via the Fouriertransform), the ω-p domain (i.e., via the Radon transform where p is avector ray parameter), the τ-p domain (i.e., the slant stack domain),etc., as those terms and transforms are known and understood in the art.Thus, although it is preferred that the 2D transform that is appliedwithin each frequency (or projection coefficient) shot or receiver slicebe a 2D Fourier transform, in some preferred embodiments the 2Dtransform could be another sort of projection transform of the sortdiscussed previously.

Additionally, it should be noted and remembered that although a typicalFourier time-direction transform of a collection of seismic tracesresults in horizontal constant frequency (or other projectioncoefficient) slices (i.e., the resulting transformed seismic traces allhave their Fourier coefficients at the same “depth” or relative arraylocation) that has been done for purposes of clarity and computationalconvenience only. Thus, when a “slice,” a “frequency slice” or a“constant frequency slice” are referred to herein those terms should bebroadly construed to include any assemblage of Fourier coefficientscorresponding to a single frequency, no matter how they happen to bearranged in memory or how they might happen to be stored on a computerreadable storage medium.

While the inventive device has been described and illustrated herein byreference to certain preferred embodiments in relation to the drawingsattached hereto, various changes and further modifications, apart fromthose shown or suggested herein, may be made therein by those skilled inthe art, without departing from the spirit of the inventive concept, thescope of which is to be determined by the following claims.

1. A method of geophysical exploration, comprising the steps of: a.selecting a plurality of seismic traces acquired over a subsurfaceexploration target, wherein each of said plurality of seismic traces hasa location associated therewith, and, wherein said plurality of seismictraces has at least one multiply reflected event therein; b. calculatinga Fourier transform of said plurality of seismic traces, therebycreating a plurality of transformed seismic traces; c. selecting anoutput trace location; d. accessing a shot record comprised of a firstplurality of said transformed seismic traces, wherein at least one ofsaid transformed seismic traces within said shot record has a locationthat is proximate to said output trace location; e. accessing a receiverrecord comprised of a second plurality of said transformed seismictraces, wherein at least one of said transformed seismic traces withinsaid receiver record has a receiver location that is proximate to saidoutput trace location; f. selecting a frequency slice within said shotrecord at least approximately corresponding to a selected frequencyvalue, thereby selecting a shot slice; g. selecting a frequency slicewithin said receiver record at least approximately corresponding to saidselected frequency value, thereby selecting a receiver slice; h.calculating a 2D Fourier transform of said shot slice, thereby obtaininga transformed shot slice; i. calculating a 2D Fourier transform of saidreceiver slice, thereby obtaining a transformed receiver slice; j. usingat least said transformed shot slice and said transformed receiver sliceto obtain an estimate of a Fourier coefficient of at least one of saidat least one multiply reflected event at said output trace location; k.performing at least steps (f) though (j) at least twice for at least twodifferent selected frequency values, thereby obtaining a plurality ofestimates of Fourier coefficients of at least one of said at least onemultiply reflected event at said output location; l. performing aninverse transform of said plurality of estimates of Fouriercoefficients, thereby obtaining an estimate of at least one of said atleast one multiply reflected event at said output location; and, m.using said estimate of at least one of said at least one multiplyreflected event at said output location to explore for hydrocarbonswithin said subsurface exploration target.
 2. The method of geophysicalexploration according to claim 1, wherein step (j) comprises the stepof: (j1) using at least said transformed shot slice and said transformedreceiver slice to obtain an estimate of a Fourier coefficient of atleast one of said at least one multiply reflected event at said outputtrace location, wherein said estimate of a Fourier coefficient isdetermined according to the equation:${{M\left( {x_{g},\left. y_{g} \middle| x_{s} \right.,{y_{s};w}} \right)} = {\sum\limits_{k_{x}}\; {\sum\limits_{k_{y}}\; {{D\left( {x_{g},\left. y_{g} \middle| k_{x} \right.,{k_{y};w}} \right)}2\; k_{z}{D\left( {k_{x},\left. k_{y} \middle| x_{s} \right.,{y_{s};w}} \right)}}}}},$where, w is said selected frequency value, M(x_(g),y_(g)|x_(s), y_(s);w) is said estimate of said Fourier coefficient at said selectedfrequency and at said output trace location, D(x_(g), y_(g)|k_(x),k_(y); w;) is a 2D Fourier transform coefficient from said transformedshot slice, D(k_(x), k_(y)|x_(s), y_(s);w) is a 2D Fourier transformcoefficient from said transformed receiver slice, 2ik_(z) is anobliquity factor, x_(g) is an X receiver coordinate, y_(g) is a Yreceiver coordinate, x_(s) is an X shot coordinate, y_(s) is a Y shotcoordinate, k_(x) is an X wavenumber, and, k_(y) is a Y wavenumber. 3.The method of geophysical exploration according to claim 1, wherein step(d) comprises the steps of: (d1) sorting said transformed seismic tracesinto a shot order, and, (d2) reading a shot record from said sortedtransformed seismic traces, wherein said read shot record is comprisedof a first plurality of said transformed seismic traces, and wherein atleast one of said transformed seismic traces within said shot record hasa location that is proximate to said output trace location.
 4. Themethod of geophysical exploration according to claim 1, wherein step (m)comprises the steps of: (m1) numerically subtracting said estimate of atleast one of said at least one multiply reflect event from at least oneof said plurality of seismic traces, thereby forming a multipleattenuated at least one seismic trace, and, (m2) using said at least onemultiple attenuated seismic trace to explore for hydrocarbons withinsaid subsurface exploration target.
 5. A method, comprising the stepsof: a. accessing a plurality of seismic traces acquired over asubsurface exploration target, at least one of said plurality of seismictraces containing a representation of at least one multiply reflectedsignal there; b. calculating a Fourier transform in time of saidplurality of seismic traces; c. accessing a shot record comprised of atleast a first portion of said plurality of Fourier transformed seismictraces, said shot record containing at least one trace proximate to aselected output trace location; d. selecting a plurality of shot sliceswithin said shot record, wherein each of said shot slices corresponds toa different Fourier frequency; e. calculating a 2D Fourier transform ofeach of said shot slices, thereby obtaining a plurality of transformedshot slices; f. accessing a receiver record comprised of at least asecond portion of said plurality of Fourier transformed seismic traces,said receiver record containing at least one trace proximate to saidselected output trace location; g. for each of said plurality oftransformed shot slices, selecting a corresponding receiver slice fromwithin said receiver record, wherein each shot slice and saidcorresponding receiver slice is associated with a same Fourierfrequency; h. calculating a 2D Fourier transform of each of saidreceiver slices, thereby obtaining a plurality of transformed receiverslices; i. using at least said plurality of transformed shot slices andsaid plurality of transformed receiver slices to obtain an estimate ofat least one of said at least one multiply reflected signal; and, j.writing to computer readable media said estimate of at least one of saidat least at least one multiply reflected signal.
 6. The method ofgeophysical exploration according to claim 5, further comprising thesteps of: k. reading at least a portion of said written estimate of atleast one of said at least one multiply reflected signal at said outputlocation, and, l. using at least said read portion of said estimate ofat least one of said at least one multiply reflected signal at saidoutput location to explore for hydrocarbons within said subsurfaceexploration target.
 7. The method of geophysical exploration accordingto claim 5, wherein step (f) comprises the steps of: (f1) sorting saidtransformed seismic traces into a receiver order, and, (f2) reading areceiver record from said sorted transformed seismic traces, saidreceiver recording being comprised of at least a second portion of saidtransformed seismic traces, and, said receiver record containing atleast one transformed seismic trace with a location proximate to saidoutput trace location.
 8. The method of geophysical explorationaccording to claim 5, wherein step (g) comprises the step of: (g1) usingat least said plurality of transformed shot slices and said plurality oftransformed receiver slices to obtain an estimate of at least one ofsaid at least one multiply reflected signal, wherein said estimate of atleast one of said at least one multiply reflected signal is determinedaccording to the equation:${{M\left( {x_{g},\left. y_{g} \middle| x_{s} \right.,{y_{s};w}} \right)} = {\sum\limits_{k_{x}}\; {\sum\limits_{k_{y}}\; {{D\left( {x_{g},\left. y_{g} \middle| k_{x} \right.,{k_{y};w}} \right)}2\; k_{z}{D\left( {k_{x},\left. k_{y} \middle| x_{s} \right.,{y_{s};w}} \right)}}}}},$where, w is a selected frequency, M(x_(g),y_(g)|x_(s), y_(s); w) is aFourier transform coefficient of said at least one of said at least onemultiply reflected signal at said selected frequency, D(x_(g),y_(g)|k_(x), k_(y); w;) is a Fourier transform coefficient from aselected transformed shot slice corresponding to said selectedfrequency, D(k_(x), k_(y)|x_(s), y_(s);w) is a Fourier transformcoefficient from said a selected transformed receiver slice at saidselected frequency, 2id_(z) is an obliquity factor, x_(g) is an Xreceiver coordinate, y_(g) is a Y receiver coordinate, x_(s) is an Xshot coordinate, y_(s) is a Y shot coordinate, k_(x) is an X wavenumber,and, k_(y) is a Y wavenumber.
 9. A method of geophysical explorationcomprising the steps of: a. accessing a plurality of seismic tracesacquired over a subsurface exploration target, and wherein at least oneof said plurality of seismic traces has at least one multiply reflectedevent therein; b. calculating a time direction discrete first projectiontransform of said seismic traces; c. accessing a shot record comprisedof a first plurality of said first projection transformed seismictraces; d. accessing a receiver record comprised of a second pluralityof said first projection transformed seismic traces; e. calculating aplurality of horizontal 2D discrete second projection transforms withinsaid shot record, thereby obtaining a plurality of transformed shotslices; f. for each of said plurality of transformed shot slices,selecting a corresponding receiver slice, said transformed shot sliceand said corresponding receiver slice having a same first projectiontransform coefficient associated therewith; g. using at least saidplurality of transformed shot slices and said corresponding plurality oftransformed receiver slices to obtain an estimate of at least one ofsaid at least one multiply reflect event; and, h. using said obtainedestimate of at least one of said at least one multiply reflected eventto explore said subsurface exploration target.
 10. The method ofgeophysical exploration according to claim 9 wherein said first discreteprojection transform is a Fourier transform and wherein said firstprojection transform coefficient is a frequency.
 11. The method ofgeophysical exploration according to claim 9, wherein said firstprojection transform and said second projection transform are a sameprojection transform.
 12. The method of geophysical explorationaccording to claim 11, wherein said first projection transform and saidsecond projection transform are a Fourier transform.
 13. The method ofgeophysical exploration according to claim 9 wherein said first discreteprojection transform is selected from a group consisting of a Fouriertransform, a fast Fourier transform, a discrete Fourier transform, aRadon transform, a cosine transform, a sine transform, a tau-ptransform, and a Walsh transform.
 14. The method of geophysicalexploration according to claim 9, wherein step (h) comprises the stepsof: (h1) numerically subtracting said estimate of at least one of saidat least one multiply reflect event from at least one seismic trace,thereby forming at least one multiple attenuated seismic trace, and,(h2) using said at least one multiple attenuated seismic trace toexplore for hydrocarbons within said subsurface exploration target. 15.A method of geophysical exploration comprising the steps of: a.accessing a plurality of seismic traces acquired over a subsurfaceexploration target, and wherein at least one of said plurality ofseismic traces has at least one multiply reflected event therein; b.calculating a time direction Fourier transform of each of said pluralityof said seismic traces; c. accessing a shot record comprised of a firstplurality of said Fourier transformed seismic traces; d. accessing areceiver record comprised of a second plurality of said Fouriertransformed seismic traces; e. calculating a plurality of horizontal 2DFourier transforms within said shot record, thereby obtaining aplurality of transformed shot slices; f. for each of said plurality oftransformed shot slices, selecting a corresponding receiver slice, saidtransformed shot slice and said corresponding receiver slice having asame frequency associated therewith; g. using at least said plurality oftransformed shot slices and said corresponding plurality of transformedreceiver slices to obtain an estimate of at least one of said at leastone multiply reflect event; and, h. using said obtained estimate of atleast one of said at least one multiply reflected event to explorewithin the subsurface exploration target.
 16. The method of geophysicalexploration according to claim 15, wherein step (g) comprises the stepsof: (g1) using at least said plurality of transformed shot slices andsaid corresponding plurality of transformed receiver slices to obtain anestimate of a corresponding plurality of Fourier transform coefficients,and, (g2) performing an inverse Fourier transform of said plurality ofFourier transform coefficients, thereby obtaining an estimate of atleast one of said at least one multiply reflected event.
 17. The methodof geophysical exploration according to claim 15, wherein step (g)comprises the steps of: (g1) selecting an output trace location, (g2)selecting a transformed shot slice at least approximately correspondingto a selected frequency value, (g3) selecting a transformed frequencyslice at least approximately corresponding to said selected frequencyvalue, (g4) using at least said selected transformed shot slice and saidselected transformed receiver slice to obtain an estimate of a Fouriercoefficient of at least one of said at least one multiply reflectedevent at said output trace location, (g5) performing at least steps (g1)though (g4) at least twice for at least two different selected frequencyvalues, thereby obtaining a plurality of estimates of Fouriercoefficients of at least one of said at least one multiply reflectedevent at said output location, and, (g6) performing an inverse transformof said plurality of estimates of Fourier coefficients, therebyobtaining an estimate of at least one of said at least one multiplyreflected event at said output location.
 18. The method of geophysicalexploration according to claim 15, wherein step (g) comprises the stepof: (g1) using at least said plurality of transformed shot slices andsaid plurality of transformed receiver slices to obtain an estimate ofat least one of said at least one multiply reflected event, wherein saidestimate of at least one of said at least one multiply reflected eventis determined according to the equation:${{M\left( {x_{g},\left. y_{g} \middle| x_{s} \right.,{y_{s};w}} \right)} = {\sum\limits_{k_{x}}\; {\sum\limits_{k_{y}}\; {{D\left( {x_{g},\left. y_{g} \middle| k_{x} \right.,{k_{y};w}} \right)}2\; k_{z}{D\left( {k_{x},\left. k_{y} \middle| x_{s} \right.,{y_{s};w}} \right)}}}}},$where, w is a selected frequency, M(x_(g),y_(g)|x_(s), y_(s); w) is aFourier transform coefficient of said at least one of said at least onemultiply reflected signal at said selected frequency, 2ik_(z) is anobliquity factor, x_(g) is an X receiver coordinate, y_(g) is a Yreceiver coordinate, x_(s) is an X shot coordinate, y_(s) is a Y shotcoordinate, k_(x) is an X wavenumber, and, k_(y) is a Y wavenumber,D(x_(g), y_(g)|k_(x),k_(y); w;) is a Fourier transform coefficient froma selected transformed shot slice corresponding to said selectedfrequency, and, D(k_(x), k_(y)|x_(s), y_(s);w) is a Fourier transformcoefficient from a selected transformed receiver slice corresponding tosaid selected frequency.